#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

// ========================= 高精度小数（模拟Java BigDecimal） =========================
// 简化版BigDecimal：用字符串存储数字，支持基础运算（可替换为GMP库提升精度）
typedef struct {
    char* value;       // 存储数字字符串（如"3.1415926"）
    int scale;         // 小数位数
    int sign;          // 符号：1正，-1负
} BigDecimal;

// 全局常量：最大允许误差（对应Java MAX_ERROR）
static const BigDecimal MAX_ERROR = {.value = "0.000001", .scale = 6, .sign = 1};

// 初始化BigDecimal（字符串转结构体）
BigDecimal bd_init(const char* str) {
    BigDecimal bd;
    bd.sign = 1;
    bd.scale = 0;
    bd.value = (char*)malloc(strlen(str) + 1);
    strcpy(bd.value, str);

    // 解析符号和小数位数
    char* dot = strchr(bd.value, '.');
    if (dot) {
        bd.scale = strlen(dot + 1);
        *dot = '\0'; // 分离整数和小数部分（简化处理）
    }
    if (bd.value[0] == '-') {
        bd.sign = -1;
        memmove(bd.value, bd.value + 1, strlen(bd.value));
    }
    return bd;
}

// 释放BigDecimal内存
void bd_free(BigDecimal bd) {
    free(bd.value);
}

// BigDecimal绝对值
BigDecimal bd_abs(BigDecimal bd) {
    bd.sign = 1;
    return bd;
}

// BigDecimal减法（a - b）
BigDecimal bd_sub(BigDecimal a, BigDecimal b) {
    // 简化实现：转换为double计算（高精度场景需替换为字符串运算）
    double a_val = atof(a.value) * a.sign;
    double b_val = atof(b.value) * b.sign;
    double res = a_val - b_val;

    char buf[64];
    snprintf(buf, sizeof(buf), "%.10f", res);
    return bd_init(buf);
}

// BigDecimal比较（a.compareTo(b)）：1(a>b)/0(a==b)/-1(a<b)
int bd_compare(BigDecimal a, BigDecimal b) {
    double a_val = atof(a.value) * a.sign;
    double b_val = atof(b.value) * b.sign;
    if (a_val > b_val) return 1;
    if (a_val < b_val) return -1;
    return 0;
}

// BigDecimal最大值（模拟Arrays.stream().max()）
BigDecimal bd_max(BigDecimal* arr, int len) {
    BigDecimal max = arr[0];
    for (int i = 1; i < len; i++) {
        if (bd_compare(arr[i], max) > 0) {
            max = arr[i];
        }
    }
    return max;
}

// ========================= 核心数据结构 =========================
// 区间结构体（对应Java Region类）
typedef struct {
    BigDecimal start;  // 区间左边界
    BigDecimal end;    // 区间右边界
} Region;

// 误差结果结构体（对应Java Error类）
typedef struct {
    BigDecimal point;  // 误差最大点
    BigDecimal error;  // 最大误差值
} ErrorPoint;

// 线性方程组结构体（对应Java EquationGroup类）
typedef struct {
    int rows;          // 行数（n+2）
    int cols;          // 列数（n+3：系数+常数项）
    double** matrix;   // 增广矩阵
} EquationGroup;

// ========================= 工具函数 =========================
// 生成区间内随机值（对应Java initRandomValue）
BigDecimal* init_random_value(int n, BigDecimal min, BigDecimal max, int* out_len) {
    *out_len = n + 2; // n+2个随机点（原逻辑隐含）
    BigDecimal* values = (BigDecimal*)malloc(*out_len * sizeof(BigDecimal));

    double min_val = atof(min.value) * min.sign;
    double max_val = atof(max.value) * max.sign;
    double range = max_val - min_val;

    srand((unsigned int)time(NULL));
    for (int i = 0; i < *out_len; i++) {
        double r = (double)rand() / RAND_MAX; // 0~1随机数
        double val = min_val + range * r;

        char buf[64];
        snprintf(buf, sizeof(buf), "%.10f", val);
        values[i] = bd_init(buf);
    }
    return values;
}

// 构建线性方程组（对应Java getEquationGroup）
EquationGroup* build_equation_group(int n, BigDecimal (*target)(BigDecimal), BigDecimal* points, int points_len) {
    if (points_len != n + 2) return NULL;

    EquationGroup* eq = (EquationGroup*)malloc(sizeof(EquationGroup));
    eq->rows = n + 2;
    eq->cols = n + 3; // 系数列(n+2) + 常数项列

    // 初始化增广矩阵
    eq->matrix = (double**)malloc(eq->rows * sizeof(double*));
    for (int i = 0; i < eq->rows; i++) {
        eq->matrix[i] = (double*)calloc(eq->cols, sizeof(double));

        // 计算目标函数值
        BigDecimal x = points[i];
        BigDecimal f_x = target(x);
        double fx_val = atof(f_x.value) * f_x.sign;

        // 构建方程：p(x) + (-1)^i * E = f(x)
        // p(x) = a0 + a1*x + a2*x² + ... + an*xⁿ
        double x_val = atof(x.value) * x.sign;
        double pow_x = 1.0; // x^0

        for (int j = 0; j <= n; j++) {
            eq->matrix[i][j] = pow_x; // 系数项：x^j
            pow_x *= x_val;
        }

        // 误差项：(-1)^i
        eq->matrix[i][n+1] = (i % 2 == 0) ? 1.0 : -1.0;
        // 常数项：f(x)
        eq->matrix[i][n+2] = fx_val;

        bd_free(f_x); // 释放目标函数返回值
    }
    return eq;
}

// 高斯消元求解线性方程组（对应Java equationGroup.resolve()）
BigDecimal* resolve_equation_group(EquationGroup* eq, int n, int* out_len) {
    *out_len = n + 2; // 解的长度：n+1个系数 + 1个误差值
    BigDecimal* res = (BigDecimal*)malloc(*out_len * sizeof(BigDecimal));

    int rows = eq->rows;
    int cols = eq->cols;
    double** mat = eq->matrix;

    // 高斯消元核心逻辑（简化版，仅处理唯一解场景）
    for (int i = 0; i < rows; i++) {
        // 找主元
        int pivot = i;
        for (int j = i; j < rows; j++) {
            if (fabs(mat[j][i]) > fabs(mat[pivot][i])) {
                pivot = j;
            }
        }
        // 交换行
        double* temp = mat[i];
        mat[i] = mat[pivot];
        mat[pivot] = temp;

        // 归一化主元行
        double div = mat[i][i];
        for (int j = i; j < cols; j++) {
            mat[i][j] /= div;
        }

        // 消去其他行
        for (int j = 0; j < rows; j++) {
            if (j != i && fabs(mat[j][i]) > 1e-10) {
                double factor = mat[j][i];
                for (int k = i; k < cols; k++) {
                    mat[j][k] -= factor * mat[i][k];
                }
            }
        }
    }

    // 提取解（系数+误差值）
    for (int i = 0; i < *out_len; i++) {
        char buf[64];
        snprintf(buf, sizeof(buf), "%.10f", mat[i][cols-1]);
        res[i] = bd_init(buf);
    }

    // 释放方程组内存
    for (int i = 0; i < rows; i++) {
        free(mat[i]);
    }
    free(mat);
    free(eq);

    return res;
}

// 生成区间分割（对应Java getRegion）
Region* get_region(BigDecimal* points, int points_len, BigDecimal min, BigDecimal max, int* out_len) {
    // 排序随机点
    double* vals = (double*)malloc(points_len * sizeof(double));
    for (int i = 0; i < points_len; i++) {
        vals[i] = atof(points[i].value) * points[i].sign;
    }
    for (int i = 0; i < points_len-1; i++) {
        for (int j = 0; j < points_len-i-1; j++) {
            if (vals[j] > vals[j+1]) {
                double temp = vals[j];
                vals[j] = vals[j+1];
                vals[j+1] = temp;
            }
        }
    }

    // 生成n+2个区间（原逻辑：10个分界点→11个区域）
    *out_len = points_len;
    Region* regions = (Region*)malloc(*out_len * sizeof(Region));

    // 第一个区间：[min, 第一个点]
    regions[0].start = min;
    regions[0].end = points[0];

    // 中间区间：[第i个点, 第i+1个点]
    for (int i = 1; i < points_len-1; i++) {
        regions[i].start = points[i-1];
        regions[i].end = points[i];
    }

    // 最后一个区间：[最后一个点, max]
    regions[points_len-1].start = points[points_len-1];
    regions[points_len-1].end = max;

    free(vals);
    return regions;
}

// 二分法找区间内最大误差点（对应Java getMaxErrorPoint）
ErrorPoint get_max_error_point(Region region, BigDecimal (*target)(BigDecimal), BigDecimal* horners, int n) {
    ErrorPoint ep;
    BigDecimal left = region.start;
    BigDecimal right = region.end;

    // 二分迭代次数
    const int iter = 50;
    for (int i = 0; i < iter; i++) {
        BigDecimal mid = bd_sub(right, left);
        // 简化：mid = (left + right)/2
        double l_val = atof(left.value) * left.sign;
        double r_val = atof(right.value) * right.sign;
        double m_val = (l_val + r_val) / 2.0;

        char buf[64];
        snprintf(buf, sizeof(buf), "%.10f", m_val);
        BigDecimal mid_bd = bd_init(buf);

        // 计算误差：target(x) - p(x)
        BigDecimal f_x = target(mid_bd);
        double fx_val = atof(f_x.value) * f_x.sign;

        // 计算多项式值 p(x) = a0 + a1*x + ... + an*xⁿ
        double px_val = 0.0;
        double pow_x = 1.0;
        double x_val = m_val;
        for (int j = 0; j <= n; j++) {
            double a_j = atof(horners[j].value) * horners[j].sign;
            px_val += a_j * pow_x;
            pow_x *= x_val;
        }

        double error = fabs(fx_val - px_val);
        snprintf(buf, sizeof(buf), "%.10f", error);
        ep.error = bd_init(buf);
        ep.point = mid_bd;

        // 二分更新（简化：向误差更大方向移动）
        // 高精度场景需对比mid左右误差
        left = mid_bd;
        bd_free(f_x);
    }

    return ep;
}

// ========================= 核心Remez算法 =========================
BigDecimal* remez(int n, BigDecimal min, BigDecimal max, BigDecimal (*target)(BigDecimal), int* out_len) {
    *out_len = n + 1; // 返回n+1个多项式系数
    BigDecimal* horners = (BigDecimal*)malloc(*out_len * sizeof(BigDecimal));

    BigDecimal step1_error = bd_init("0");
    BigDecimal step2_error = bd_init("0");
    int first_iter = 1;

    // 迭代优化：误差差的绝对值 > MAX_ERROR
    while (first_iter || bd_compare(bd_abs(bd_sub(step1_error, step2_error)), MAX_ERROR) > 0) {
        first_iter = 0;

        // Step 1：生成随机点，构建并求解方程组
        int rand_len = 0;
        BigDecimal* random_values = init_random_value(n, min, max, &rand_len);

        EquationGroup* eq = build_equation_group(n, target, random_values, rand_len);
        int resolve_len = 0;
        BigDecimal* resolve = resolve_equation_group(eq, n, &resolve_len);

        // 第一步最大误差（解的最后一个值）
        bd_free(step1_error);
        step1_error = resolve[resolve_len - 1];

        // Step 2：生成区间，找每个区间最大误差点
        int region_len = 0;
        Region* regions = get_region(random_values, min, max, &region_len);

        // 复制多项式系数到horners
        for (int i = 0; i <= n; i++) {
            horners[i] = resolve[i];
        }

        // 计算每个区间的最大误差
        BigDecimal* errors = (BigDecimal*)malloc(region_len * sizeof(BigDecimal));
        for (int i = 0; i < region_len; i++) {
            ErrorPoint ep = get_max_error_point(regions[i], target, horners, n);
            random_values[i] = ep.point;
            errors[i] = ep.error;
        }

        // 第二步最大误差
        bd_free(step2_error);
        step2_error = bd_max(errors, region_len);

        // 释放临时内存
        for (int i = 0; i < rand_len; i++) {
            bd_free(random_values[i]);
        }
        free(random_values);
        for (int i = 0; i < resolve_len - 1; i++) { // 保留最后一个误差值
            bd_free(resolve[i]);
        }
        free(resolve);
        for (int i = 0; i < region_len; i++) {
            bd_free(regions[i].start);
            bd_free(regions[i].end);
            bd_free(errors[i]);
        }
        free(regions);
        free(errors);
    }

    // 释放临时误差变量
    bd_free(step1_error);
    bd_free(step2_error);

    return horners;
}

// ========================= 测试示例（目标函数：sin(x)） =========================
BigDecimal target_sin(BigDecimal x) {
    double x_val = atof(x.value) * x.sign;
    double sin_val = sin(x_val);

    char buf[64];
    snprintf(buf, sizeof(buf), "%.10f", sin_val);
    return bd_init(buf);
}

int main() {
    // 测试参数：3次多项式，区间[0, π]，目标函数sin(x)
    int n = 3;
    BigDecimal min = bd_init("0.0");
    BigDecimal max = bd_init("3.1415926535");

    int out_len = 0;
    BigDecimal* horners = remez(n, min, max, target_sin, &out_len);

    // 打印多项式系数
    printf("Remez算法优化后的%d次多项式系数：\n", n);
    for (int i = 0; i < out_len; i++) {
        printf("a%d = %s\n", i, horners[i].value);
        bd_free(horners[i]);
    }
    free(horners);

    bd_free(min);
    bd_free(max);
    return 0;
}